\documentclass[12pt, a4paper]{article}
\usepackage{geometry}
\geometry{
	a4paper,
	left=12.7 mm,
	right=12.7 mm,
	top=12.7 mm,
	bottom=12.7 mm,
}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{physics}
\graphicspath{ {./} }
\usepackage{ctex}

\newcommand{\bvec}[1]{\ensuremath{\mathbf{#1}}}
\newcommand{\formula}[1]{\text{式} \ref{#1} }
\newcommand{\matrx}[1]
{
	\ensuremath
	{
		\left (
		\begin{matrix}
			#1
		\end{matrix}
		\right)
	}
}

\begin{document}
	\setcounter{section}{-1}
	
	\section{让我们问问神奇AI}
	
	“
	好的，我来给你讲一讲PDE（偏微分方程）的有限差分法。
	
	想象一下，你有一个很大的画纸，上面画了很多很多的小格子。每个小格子上都有一个数字，这个数字代表了这个小格子的温度。现在，我们想要知道这些小格子的温度是如何变化的。
	
	有限差分法就像是一个游戏规则，它可以帮助我们预测这些小格子的温度变化。
	这个游戏规则是这样的：我们看看每个小格子和它周围的小格子（比如上面、下面、左边、右边的小格子），然后根据这些小格子的温度来猜测中间的小格子的温度会如何变化。
	
	具体来说，我们可以用一个简单的公式来表示这个猜测。比如说，如果我们想知道中间的小格子的温度会如何上升或下降，我们可以看看它和上下左右格子的温度差。
	然后，我们可以用这些温度差来猜测中间的小格子的温度变化。
	
	我们把这张画满小格子的纸放在计算机里，计算机就会用这个游戏规则来帮助我们知道每个小格子的温度变化。
	计算机可以非常非常快地做这个游戏，所以它可以告诉我们很多很多小格子的温度变化。这样，我们就可以知道整个画纸上的温度是如何变化的了。
	
	这个方法就是有限差分法，它是用计算机来帮助我们解决一个非常复杂的数学问题的一种简单而有效的方法。
	”
	
	---摘录自AI回复
	
	\section{PDE方程}
	\footnote{本文是Sauer《数值分析》的学习笔记} 
	ODE方程中，函数只有一个自变量，因此ODE方程只包括一个自变量的导数。
	$$y=f(x) ~ \text{such that} ~ \dv{y}{t} = f(y,t)$$
	而PDE方程中，函数有多个自变量，因而PDE方程包括函数对各个自变量的各阶偏导数。
	比如我们熟知的扩散方程是一个PDE，其包含$u$（代表温度或者浓度等物理量）关于空间与时间的导数：
	\begin{equation}
		u=u(x,y,z,t) ~ \text{such that} ~ \pdv{u}{t} = D \laplacian{u}
	\end{equation}
	其中 $\laplacian{u} = \pdv[2]{u}{x}+\pdv[2]{u}{y}+\pdv[2]{u}{z}$。
	
	不少物理方程都是PDE；同时，PDE解起来\textsl{也}比ODE难多了。
	\textsl{为什么重点、考点和难点总是出现在一起！}
	
	\newpage
	\section{显式法}
	
	\subsection{扩散方程}
	
	\begin{figure}[h]
		\centering
		\includegraphics[width=0.5\linewidth]{diffiter}
		\caption{扩散方程迭代示意图}
		\label{fig:diffiter}
	\end{figure}
	
	我们以一维扩散方程为例，探讨如何求解PDE。
	\begin{equation}
		u = u(x, t), 
		\left \{
		\begin{aligned}
			\pdv{u}{t} &= D \pdv[2]{u}{x}\qquad \text{扩散方程}\\
			u(x,t_0) &= f(x) \qquad \text{初始条件}\\
			u(x_0,t) &= l(t) \qquad \text{左侧边界条件}\\
			u(x_n,t) &= r(t) \qquad \text{右侧边界条件}\\
			x & \in [x_0, x_n]\\
			t & \ge t_0 \\
		\end{aligned}
		\right.
	\end{equation}
	和ODE一样，我们得分别离散化时间和空间：
	我们取离散时刻为$t_k = t_0 + k \Delta t, k=0,1,2,3,...$, 离散空间为$x_i = x_0 + i \Delta x, i = 0,1,2,3,..., n$；
	随后，分别求解在这些离散时刻、空间上的场$u = u(x_i, t_k)$。
	为简化书写，我们记$u^{(k)}_{i} = u(x_i, t_k)$。	
	此外，我们还需要知道初始条件与边界条件。
	
	仿照在ODE中的操作，我们分别离散化两个导数项：
	\begin{equation}
		\begin{aligned}
			\pdv{u(x,t)}{t} &= \frac{u(x,t+\Delta t) - u(x,t)}{\Delta t} \\
			\pdv[2]{u(x, t)}{x} &= \frac{u(x+\Delta x, t) - 2u(x, t) + u(x-\Delta x, t)}{(\Delta x)^2} \\
		\end{aligned}
	\end{equation}
	那么扩散方程变为离散形式：
	\begin{equation}
		\begin{aligned}
			& \pdv{u}{t} = D \pdv[2]{u}{x} \\
			\Rightarrow & \frac{u^{(k+1)}_i - u^{(k)}_i }{\Delta t} = D \frac{u^{(k)}_{i+1} - 2u^{(k)}_{i} + u^{(k)}_{i-1}}{(\Delta x)^2}\\
			\Rightarrow & u^{(k+1)}_{i}  = u^{(k)}_{i} + \frac{D \Delta t}{(\Delta x)^2} [u^{(k)}_{i+1}  - 2u^{(k)}_{i}  +u^{(k)}_{i-1} ]\\
		\end{aligned}
	\end{equation}
	这个公式适用于 $i = 1,2,3,...,n-1$；
	$u_0$与$u_n$将直接由边界条件确定。
	总而言之，
	\begin{equation} \label{eq_numdif}
		\boxed
		{
			\begin{aligned}
				u^{(k+1)}_{i}  &= u^{(k)}_{i} + \frac{D \Delta t}{(\Delta x)^2} [u^{(k)}_{i+1}  - 2u^{(k)}_{i}  +u^{(k)}_{i-1} ] \qquad i=1,2,3,...,n-1 \\
				u^{(k)}_0 &= l(t_k), u^{(k)}_n = r(t_k) \\
			\end{aligned}
			\qquad k = 0,1,2,3,..., 
		}
	\end{equation}
	
	\newpage
	\subsection{波动方程}
	
	\begin{figure}[h]
		\centering
		\includegraphics[width=0.5\linewidth]{waveiter}
		\caption{波动方程迭代示意图}
		\label{fig:waveiter}
	\end{figure}
	
	随后我们处理一维波动方程，方法也是大同小异：
	\begin{equation}
		u = u(x, t), 
		\left \{
		\begin{aligned}
			\pdv[2]{u}{t} &= c^2 \pdv[2]{u}{x}\qquad \text{波动方程}\\
			u(x,t_0) &= f(x) \qquad \text{初始条件}\\
			\pdv{u}{t}|_{t=t_0} &= g(x) \qquad \text{速度初始条件}\\
			u(x_0,t) &= l(t) \qquad \text{左侧边界条件}\\
			u(x_n,t) &= r(t) \qquad \text{右侧边界条件}\\
			x & \in (x_0, x_n)\\
			t & \ge t_0 \\
		\end{aligned}
		\right.
	\end{equation}
	现在，$u$关于时间和空间的导数都是二阶的。
	和二阶ODE一样，我们不仅要知道初始条件，还得知道初始速度条件。
	我们仍然分别离散化二阶导数：
	\begin{equation}
		\begin{aligned}
			\pdv[2]{u(x,t)}{t} &= \frac{u(x, t + \Delta t) - 2u(x, t) + u(x, t- \Delta t)}{(\Delta t)^2} \\
			\pdv[2]{u(x, t)}{x} &= \frac{u(x+\Delta x, t) - 2u(x, t) + u(x-\Delta x, t)}{(\Delta x)^2} \\
		\end{aligned}
	\end{equation}
	因此，
	\begin{equation}
		\begin{aligned}
			& \pdv[2]{u}{t} = c^2 \pdv[2]{u}{x} \\
			\Rightarrow & \frac{u^{(k+1)}_{i} - 2u^{(k)}_{i} + u^{(k-1)}_{i}}{(\Delta t)^2} = c^2 \frac{u^{(k)}_{i+1} - 2u^{(k)}_{i} + u^{(k)}_{i-1}}{(\Delta x)^2} \\
			\Rightarrow & u^{(k+1)}_i = 2u^{(k)}_i - u^{(k-1)}_i+ \left( \frac{c \Delta t}{\Delta x} \right)^2 (u^{(k)}_{i+1} - 2u^{(k)}_{i} + u^{(k)}_{i-1})\\
		\end{aligned}
	\end{equation}
	当我们取$k=0$计算$u^{(1)}$时，我们会遇到烦人的\textsl{Missing No.} $u^{(-1)}$。
	解决方法参考ODE的附件，即我们要使用不同的方法计算$u^{(1)}$。
	\begin{equation}  \label{eq_numwav}
		\boxed
		{
			\begin{aligned}
				u^{(k+1)}_i &= 2u^{(k)}_i - u^{(k-1)}_i+ \frac{c^2 (\Delta t)^2 }{(\Delta x)^2} (u^{(k)}_{i+1} - 2u^{(k)}_{i} + u^{(k)}_{i-1}) \qquad k = 1,2,3,..; i=1,2,3,...,n-1\\
				u^{(1)}_i & = 1/2 \left(2u^{(0)}_i + 2g(x_i)\Delta t+  \left( \frac{c \Delta t}{\Delta x} \right)^2 (u^{(0)}_{i+1} - 2u^{(0)}_{i} + u^{(0)}_{i-1} ) \right) \qquad i=1,2,3,...,n-1 \\
				u^{(k)}_0 &= l(t_k), u^{(k)}_n = r(t_k) \\
			\end{aligned}
		}
	\end{equation}	
	
	\newpage
	\section{边界}
	
	\subsection*{第一类边界}
	
	上文中，我们直接设置了$u$在边界处的值，这种边界称为\textbf{第一类边界}：
	\begin{equation}
		u(x_0,t) = l(t) \qquad u(x_n,t) = l(t)
	\end{equation}
	
	\subsection*{第二类边界}
	我们实则还可以给定边界的导数值，称为\textbf{第二类边界}：
	\begin{equation}
		\pdv{u(x,t)}{x}|_{x=x_0} = q_l(t) \qquad \pdv{u(x,t)}{x}|_{x=x_n} = q_r(t)
	\end{equation}	
	我们使用二阶导的单侧导数形式离散化导数：
	\begin{equation}
		\pdv{u(x,t)}{x} |_{x=x_0} = \frac{- 3u(x_0,t)+ 4u( x_0 + \Delta x,t) -u(x_0 + 2\Delta x,t)}{2 \Delta x}
	\end{equation}
	整理上述等式，我们发现应该如此设置$u$的边界值以应用第二类边界：
	\begin{equation}
		u(x_0,t) = \frac{ - q_l(t)(2 \Delta x) +4u( x_0 + \Delta x,t) -u(x_0 + 2\Delta x,t)}{3}
	\end{equation}

	\subsection*{周期性边界}
	
	\textbf{周期性边界}假设左侧边界与右侧边界直接相连。
	这种假设在许多物理问题中非常常见，
	因为它能够有效减少边界效应。
	\textsl{或许地球也是一种周期性边界，你可以绕着地球转一圈而回到原点。}
	
	我们以左侧边界为例，
	在标准情况下，二阶导数的中心差分公式为：
	$$
	\pdv[2]{u(x, t)}{x} \bigg|_{x=x_0} = \frac{u(x_0+\Delta x, t) - 2u(x_0, t) + u(x_0-\Delta x, t)}{(\Delta x)^2}.
	$$
	很显然，$x_0-\Delta x$ 超出了定义域，因此 $u(x_0-\Delta x, t)$ 并不存在。
	然而，依据周期性边界的条件，我们可以假设 $u(x_0-\Delta x, t)$ 等价于右侧边界的值 $u(x_n, t)$，其中 $x_n$ 是右侧边界的坐标。
	那么我们得到：
	$$
	\pdv[2]{u(x, t)}{x} \bigg|_{x=x_0} = \frac{u(x_0+\Delta x, t) - 2u(x_0, t) + u(x_n, t)}{(\Delta x)^2}.
	$$
	类似地，对于右侧边界 $x_n$，可以将超出定义域的点 $x_n+\Delta x$ 映射回左侧边界 $x_0$。
	在编程实现中，可以利用数组滚动操作（如 NumPy 的 `np.roll` 或 PyTorch 的 `torch.roll` 函数）来高效处理周期性边界条件。
	这些函数会自动将数组的一端连接到另一端，从而避免手动处理边界映射的复杂性。

\end{document}
